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ABSTRACT 

We investigate the non-resonant, 3-D (spatial) model of the hierarchical system composed of 
point-mass stellar (or sub-stellar) binary and a low-mass companion (a circumbinary planet 
or a brown dwarf). We take into account the leading relativistic corrections to the Newtonian 
gravity. The secular model of the system relies on the expansion of the perturbing Hamiltonian 
in terms of the ratio of semi-major axes a, averaged over the mean anomalies. We found 
that the low-mass object in a distant orbit may excite large eccentricity of the inner binary 
when the mutual inclination of the orbits is larger than about of 60 deg. This is related to 
strong instability caused by a phenomenon which acts similarly to the Lidov-Kozai resonance 
(LKR). The secular system may be strongly chaotic and its dynamics unpredictable over the 
long-term time scale. Our study shows that in the Jupiter- or brown dwarf- mass regime 
of the low-massive companion, the restricted model does not properly describe the long-term 
dynamics in the vicinity of the LKR. The relativistic correction is significant for the parametric 
structure of a few families of stationary solutions in this problem, in particular, for the direct 
orbits configurations (with the mutual inclination less than 90 degrees). We found that the 
dynamics of hierarchical systems with small a ~ 0.01 may be qualitatively different in the 
realm of the Newtonian (classic) and relativistic models. This holds true even for relatively 
large masses of the secondaries. 

Key words: celestial mechanics - secular dynamics - analytical methods - stationary solu- 
tions - extrasolar planetary systems 



1 INTRODUCTION 

The extrasolar planets are discovered routinel)^] Recently, about 
of 500 low-mass companions to stars of different spectral types are 
known. Most of them are bounded to single stars. Moreover, there is 
also a growing number of planetary candidates in binaries as well 
as multi-stellar systems (see [Eggenberger 2010 for the statistical 
properties of planets in binaries). Generally, following nomencla- 
ture in (Rabl & Dvorak 1988), we can consider two classes of such 
multiple systems. In the satellite case or the S-type configuration, 
a planet revolves around one of the primaries in the binary, and the 
second primary is much more distant. In the cometary or circumbi- 
nary configuration (C-type from hereafter), the planet has a wide 
orbit around the inner, massive binary. 

Current theories of planet formation in multiple stellar sys- 
tems (e.g., Takeda et al. 2009 and references therein) show that the 
inclination of a planetary orbit to the orbital plane of the binary may 
be non-zero in the C-type and the S-type systems. It is well known 
that in the S-type configurations, when the mutual inclination of cir- 
cular orbits is larger than the critical value / crit ~ 40°, the inner orbit 
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undergoes large-amplitude oscillations of the eccentricity, which 
is in anti-phase with the mutual inclination. This dynamical phe- 
nomenon is well known as the Kozai (or Lidov-Kozai) resonance 
( [Lidov| 1962| |Kozai| 1962) . We will call it the LKR from hereafter. 
Keeping in mind the two types of possible orbital configurations, 
this instance of the LKR may be also understood as the inner LKR 
(see, e.g., Krasinskii 1 1 972| [T974] > . Actually, many authors explain 
large eccentricities of some planetary candidates due to forcing by 
a distant star or massive companion (see, e.g., |Tamuz, P., et al.,| 
2008 ; Fabrycky & Tremaine 2007). In fact, the amplification of ec- 
centricity and inclination may also appear in the C-type systems. 
The critical inclination is then ~ 60°, and it may be attributed to 
the outer LKR (see, e.g., |Krasinskii|[l972| |1974| |Migaszewski &| 
Gozdziewski 2010), and this will be addressed further in this work. 

In the literature, the problem is most often considered as the 
restricted one, which means that the planet is a mass-less particle 
not perturbing the motion of the binary. It has been studied in many 
papers, with different analytical and numerical techniques. The re- 
stricted model help us to simplify the analysis, nevertheless, the 
assumption of negligible influence of the planet on the motion of 
primaries may be not valid if the planet is large. In fact, low-mass 
objects in a few Jupiter-mass range are quite common. If the mutual 
interactions are significant, as we will show further in this work, the 
binary orbit may be strongly perturbed by a distant, relatively mas- 
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sive planet or a brown dwarf moving in inclined, wide orbit, even if 
the mass of inner companion is ten times larger than the perturber. 

In this work, we focus on the unrestricted C-type problem by 
means of the secular model in terms of the semi-major axes ratio, 
a. We assume that a is small (< 0.1), hence we focus on hierar- 
chical systems. In such a case, we face very different times scales 
of the orbital evolution. Typically, the inner binary has the orbital 
period counted in months or years but the period of the outer planet 
is at least ~ 100 times longer. Hence the short-term mean motion 
resonances (MMRs) are not present in the system. The secular evo- 
lution of the mean orbits, depending on the mutual interactions, is 
even much longer, and spans Myr time-scale. To follow the orbital 
evolution in all these time-scales in details, one could integrate the 
equations of motion numerically. Unfortunately, this direct, brute- 
force approach requires too large CPU overhead. 

Because we are primarily interested in the long-term evolution 
of the hierarchical systems, the problem may be simplified with 
the help of the averaging principle. The Hamiltonian of the hier- 
archical system may be splitted onto integrable Keplerian part of 
the inner binary, and for the perturbation part of the mutual inter- 
action with the low-mass companion. Because the system is non- 
resonant, the perturbing Hamiltonian may be averaged out over the 
mean anomalies or the mean longitudes, which play the role of the 
fast angles. After the averaging, we obtain the secular Hamiltonian 
describing the long-term evolution of the mean, slowly varying or- 
bits. To obtain the secular model, we extend a simple averaging 
scheme (the so called mixed anomalies method) in our earlier pa- 
per (Migaszewski & Gozdziewski 2008 ) to the non-coplanar case. 
The perturbing Hamiltonian is expanded in power series with re- 
spect to a, and then these series may be averaged out term by term. 
We derived such averaged expansion of the secular Hamiltonian to 
the 10-th order in a. It is a more general version of the third-order 
(octupole-level) theory, studied in the planetary context in ( |Ford| 
|et al.|2000||Lee & Pe ale 2003 ) and of the integrable, second order 
(quadrupole-level) approximation in many papers (e.g., |Harring- 



[tonTr9 68, Krasinskii |T972l|Lidov & Ziglin|1976||Ferrer &O sacar 
1994, Farago & Laskar 2010[ and references therein). 

Our work is closely related to remarkable study by 
Michtchenko et al. (2006}, to which we will refer many times, as 
well as to our earlier paper ( Migaszewski & Gozdziewski 2009a ) 
devoted to the analysis of equilibria in the three-dimensional prob- 
lem of two planets. Moreover, this work extends these papers in 
two important aspects. The averaging of the secular problem is 
done analytically, which simplifies and optimizes the computa- 
tions, helping us to avoid numerical artifacts. Here, we also con- 
sider a generalization of the Newtonian model, by accounting for 
the leading non-Newtonian point-mass corrections to the perturb- 
ing Hamiltonian, i.e., the relativistic, post-Newtonian (PN) correc- 
tion. It can be also easily averaged out over the mean anomalies. 
The PN corrections are particularly important if the frequencies 
of the slow angles, which they induce, are comparable with the 
frequencies causes by the Newtonian interactions (e.g., |Adams &| 
Laughlin 2006, Fabrycky & Tremaine 2007 ). We shown already 
(Migaszewski & Gozdziewski 2009b) that accounting for the rela- 
tivistic corrections in the co-planar case of two-planet system may 
lead to qualitatively different global view of the phase space in both 
models. Hence, the PN correction might be also important in the 3D 
problem. Indeed, as we will show below, this apparently subtle ef- 
fect induces global, qualitative changes of the structure of the phase 
space. 

It should be emphasised here, that a large number of physi- 
cal and orbital parameters fully characterising planetary configura- 



tions contradicts our desire to study the problem in possibly qual- 
itative way, with the help of particular geometric tools. Hence, we 
restrict the work to specific ranges of these parameters, focusing 
on a "typical" binary with relatively large mass ratio of the pri- 
maries, as well as the circumbinary object in Jupiter/brown-dwarf 
mass range. Moreover, considering corrections to the Newtonian 
point-mass gravity, we only consider the relativistic effects, which, 
in turn, limit the orbital parameters of the binary. The conserva- 
tive and dissipative tidal distortions are neglected here, though they 
might dominate in compact binaries, or in configurations with very 



hot- Jupiter or super-Earth planets (Fabrycky & Tremaine 
Mardli ngl [20071 |Batygin et aLl|2009| |Ragozzine & Wolf 



2007 



2009 



Mardling 2010). In the range of semi-major axes ~ 0.025 au, the 
planetary tidal bulge raises apsidal rotation of the inner orbit which 
may may reach a few degrees per year, exceeding the effects of gen- 
eral relativity and the rotational stellar quadrupole by more than an 
order of magnitude (Ragozzine & Wolf 2009 ). However, in gen- 
eral, as we explain below, the rotational distortions introduce extra- 
degrees of freedom to the model (assuming that the stellar and plan- 
etary spins may be arbitrary), that cannot be treated in terms of ge- 
ometric tools natural to investigate two-degrees of freedom Hamil- 
tonian dynamics. Still, although the tidal effects could be basically 
treated in this formalism too, it would introduce new parameters 
(bodies' radii, Love numbers), hence we postpone investigations of 
this more general and complex model in future papers. Overall, as 
we show below, in the parameter ranges investigated here (semi- 
major axis of the binary ~ 0.1-0.2 au and larger), the general rel- 
ativity is dominant over the rotational and tidal corrections to the 
mutual Newtonian interactions. Yet we shall also demonstrate that 
our results may be quite easily scaled down to the regime of masses 
and semi-major axes typical for multi-planet configurations, and in- 
vestigated in the coplanar case mostly. 

A plan of this paper is as follows. In Sections 2 and 3, we 
derive the 3D secular model of the planetary system, in that the 
mean motion resonances of low order are absent, following the co- 
planar case considered in (Migaszewski & Gozdziewski 2008). We 
try to keep the presentation self-contained, therefore we recall basic 
facts regarding the dynamics of the secular model, which might be 
found in other papers already published. Section 3 describes a test 
of the secular approximation, and recalls the notion of the so called 
representative planes of initial conditions, as well as a scheme of 
investigating families of stationary solutions in the secular model. 
Section 4 is for the analysis of the eccentricity evolution and chaotic 
dynamics. Section 5 is devoted to a parametric survey of the equi- 
libria in the classic (point-mass Newtonian) model. In Sect. 6, we 
study influence of the PN corrections on these solutions. In Con- 
clusions, we summarize the results and sketch perspectives of the 
further research. 



2 THE SECULAR 3D MODEL OF TV-BODIES 

We consider the general, spatial model of the Af-planet sys- 
tem around a central star. It may be described in terms of the 
Hamiltonian written with respect to canonical Poincare variables 
( [Michtchenko et al.|2006) in the form of = ^ ep i + ^ ert , where 



Ha 



:epl " 



N 

S V2P* 



(i) 



stands for the integrable part comprising of the direct sum of the rel- 
ative, Keplerian motions of point-mass secondaries m^ i = 1 , . . . , N, 



© 2008 RAS, MNRAS 000,[l][T8] 



The relativistic dynamics of circumbinary planets 3 



with respect to the primary mass mq. We also define the mass pa- 
rameters /u* = k 2 (mo + mi), where k is the Gauss gravitational con- 
stant, and p* = (1/mj + 1 /mo) -1 are the so called reduced masses. 
The term ^ ert stands for the perturbing function of the Keplerian 
motions. We assume that the perturbation is a sum of two terms: 



■^pert — -^G + ^GR, 



(2) 



where is related to a small Newtonian mutual interactions be- 
tween nti and my, and we assume that ^ e rt/-^ke P i *C 1. That may 
be accomplished either by keeping m; small (then we have the 
planetary regime) or permitting that secondary masses are rela- 
tively large (even comparable with the central object) and simul- 
taneously requiring large separations between particular orbits (the 
stellar regime). The term of is for the leading general relativity 
(PN) corrections to the potential of the central star and the inner- 
most companion. Here we focus on the non-resonant systems, with 
well separated orbits, hence we may neglect the relativistic post- 
Newtonian perturbations of the outer bodies. If the semi-major axes 
ratios GC;j = at/aj < 0.1 are small, the relativistic corrections for 
the distant objects are by orders of magnitude smaller than the PN 
perturbation of the inner binary. 

Following the notion of the Poincare coordinates, the per- 
turbation may be written as follows: 

k 2 m;m j 



N-l N 

i=l j>i 



VrVj 

mo 

indirect part 



(3) 



direct part 

where r ? , i= 1, ...,Af are for the position vectors of m* relative to 
the central body, p ; are for their conjugate momenta relative to the 
bary center of the whole (N-\- l)-body system, A/j = ||rj — ry|| de- 
note the relative distance between bodies / and j. 

After ( |Richardson & Kelly| [l988), or developing the PN 
Hamiltonian from the general Lagrangian in ( Brumberg 2007 J, the 
post-Newtonian potential in the PN formulation, ^Hq R ee P*^ R , 
where ^ R has the following form: 

^2 1 



^ R = YiP 4 +Y2y 
with coefficients 7i , 72 , 73 , 74 : 

(1-3V) ^ a/*(3- 



-73 



7i 



8c 2 ' 72 ' 2c 2 

where c stands for the speed of light in a vacuum, jf = k 2 (mo + 
m\ ), V = mom j / (mo + m\ ) 2 , P is the astrocentric momentum of the 
innermost secondary (normalized through the reduced mass): 



(r-P) 



73 : 



-74- 



2c 2 



74 : 



(4) 



(TV 

~2?' 



1 



4yi(v- v)v- 



2Y2 2ft 
v + ~ 



(5) 



and v = r stands for the astrocentric velocity of the innermost ob- 
ject (still, assuming that the relativistic corrections from the other 
bodies in the system are neglected). Hence, P = v with the accuracy 
of 0(c -2 ) and then the relativistic Hamiltonian is conserved up to 
the order of 0(c~ 4 ). 

It is well known that the equations of motion of the Af-body 
system with N ^ 3 are not integrable. However, making use of 
the assumptions above, we may apply the averaging proposition 
(Arnold et al. 1993 ) to remove the short order perturbations, and to 
derive the equations of the secular dynamics, governing the long- 
term evolution of the mean orbital elements. 

To perform the averaging, the perturbation must be expressed 
in terms of the canonical action-angle variables (!,(])): 



#(i,4>) = ^(i) + ^(i^), 



and we assume that ^ ert (I,(|)) ~ 8^4 epl (I), where e <C 1 is a small 
parameter. Here, we use the mass-weighted Delaunay elements 
(e.g., |Murray & Dermott|2000| ): 



U EE M U Li = p* V^T^' 

gi = CD;, Gi = Li ^1 - ef 
hi = H z , Hi = Gi cos Ii, 



(7) 



where Mi are the mean anomalies, a; stand for canonical semi- 
major axes, ei are the eccentricities, // denote inclinations, co ; are 
the arguments of pericenter, and £2; denote the longitudes of as- 
cending node. The Hamiltonian of the Af-planet system written in 
terms of the these Delaunay variables (Eq.[7} has the form of: 



N (u*) 2 (R*) 3 

• I „2 +Hert {L U l U GuguH U hi) . (8) 
i=l,...,N 



In this formulation, // are the fast angles. They may be eliminated 
through the averaging that is accomplished with: 



I r2% r2% 



.dM N . 



(9) 



i=l,. ..JV 



We should remember here that ^ ec is valid only if (1) i^ ert ~ ei^kepi 
(where 8<1 means a small parameter), and the averaging of the 
unperturbed Keplerian orbits is equivalent to performing the first 
step of the perturbation calculus ( Ferraz-Mello 2007 ), (2) there is 
no mixed resonances between the inner binary and the outer com- 
panion (e.g., between slow frequencies of the inner orbit and the 
mean motion of the outer orbit). We checked that planetary systems 
studied in this paper obey these assumptions within respective pa- 
rameter ranges. These calculations rely on the averaged model, and 
will be given below (see the end of Sect. 3). Because the secular 
Hamiltonian i^ ec does not depend on mean anomalies Mi, the con- 
jugate momenta Li are integrals of the secular problem. Obviously, 
the mean semi-major axes are also constant, hence we get rid of 
Af-degrees of freedom (DOF). 



3 AVERAGING THE 3D MODEL OF TV-BODIES 



In ([Migaszews ki & Goz dziewski 2008), we describe a simple 
scheme of the averaging the perturbing function ^ ert (Eq[2]) in 
coplanar case, which makes use of the very basic properties of 
the Keplerian motion, the mixed anomalies algorithm. This method 
may be easily applied to the 3D problem. At first, we consider the 
direct part of the mutual interaction between the planets, (Eq. 
[3}. The indirect part averages out to a constant and does not con- 
tribute to the secular dynamics (Brouwer & Clemence 1961 ). 

The secular Hamiltonian may be written as a sum of terms 
representing mutual interactions between all pairs of bodies i < j, 
where 0/ < ay. 



N-l N 
(Hie) = Y<H 



(10) 



For a particular pair of planets, we calculate the following integral: 



(6) 



Hence, the problem may be reduced to averaging the inverse of 
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the distance A/j between two particular planets over their mean 
anomalies: 



(12) 



where r* and rj must be expressed in a common reference frame 
jF. The same vectors, written in the orbital reference frames ft of 
each planet, are r/|^- = [x;,y;,0] r , and expressed in the common 

reference frame, they have the form of: = A/r ; I ^ . Here, the rota- 

i y-i 

tion matrix A/ = A/(co;,£^,/;) is the matrix product of elementary 
Eulerian rotations (Murray & Dermott 2000): 

Formulae 1 1 2 1 may be rewritten as follows: 



m 1-2- 



(13) 



Following ( Migaszewski & Gozdziewski 2008), we express the ra- 
dius vector r; of the inner body in each planetary pair with respect 
to the eccentric anomaly. The radius vector of the outer body in the 
pair is parametrised by the true anomaly. This choice implies that 
A^ 1 expanded in Taylor series with respect to a is a trigonometric 
polynomial. To compute the integral in Eq.[TT] we also change the 
integration variables: 



dMi = IfdEi, and dMj = Jjdfj, 



where auxiliary functions appear: 



Ii = l-eiCOsEt, Jj 



(l-^-) 3/2 
(l + ejco&fj) 



(14) 



Finally, the averaged mutual perturbation has the same general 
form as in the coplanar case (Migaszewski & Gozdziewski 2008): 



k 2 L 



(15) 



1=2 



although explicit expressions for are obviously different in 

the 3D model. The zeroth-order term in Eq. [15] is reduced to a 
constant and does not influence the secular equations of motion. 
Also the first order term vanishes identically. The remaining terms 
have rather complex form. In the simplest three-body sys- 
tem (i = 1,7 = 2), we may express them in the Laplace reference 
frame. In this case, AO. = ±n and G\ sin/i = G 2 sin/2 (see e.g., 
Michtchenko et al. 2006). It is also natural to introduce the mutual 
inclination, / mut = I\ + I 2 . Then the $Q -terms of the order of 2 
and 3 may be identified with the quadrupole and octupole terms, 
respectively (see, e.g., Ford et al. 2000, Lee & Peale 2003 1 Farago 
|& Las kar 2010). The quadrupole-level term is the following: 



2 + 3e- 



!5 2 

— eiD 2 cos 2coi, 
16 1 



(16) 



where Q = cos/ mut , and D x = (3Cf - l)/2, D 2 = Cj -I. The 



third order (octupoloe-level) term reads as follows: 

1,2) _ 15 , 



^D 6 eie 2 cosAG5 + 4^ 

256 3 D2e i e2COS ( 30)1 _c ° 2 ) 
525 
"512 
15 



D4D2qe2COs(3a>i +CO2) 



(17) 



+ T28 D5 ^^ 2 cos ^ c ° 1 + 0)2 ) 



45 
512' 



+ ^ D 5^1^2COS(C0 1 +C0 2 ) 



where coefficients Dj are: 



Z>3 = (l + C/)/2, D 4 = l-C h 
D 5 = -15C|+5C^ + 11C/-1, 
D 6 = (15C? + 5C?-llC/-l)/8. 



Equations [16] and[l7] are written similarly to terms appearing in the 
coplanar model [see equations (22) and (23) in (Migaszewski & 
|Gozdziewski|2008| )]. Clearly, if z mut = then D x = D 3 = D 6 = 1, 



D 2 =D 4 =D 5 = 0, and formulae ( K} 3 J) coincide with those 

ones of in the coplanar problem. 

An explicit expansion of ^ ec shows that the quadrupole-order 
term in a introduces the evolution of eccentricity e\, and in this ap- 
proximation, the outer eccentric e 2 is constant. The variation of the 
outer eccentricity may be only introduced through the third order 
(octupole) and higher terms. Indeed, up to the quadrupole approxi- 
mation, the secular Hamiltonian does not depend on 002 (the cyclic 
angle), and the eccentricity of the outer body becomes an additional 
integral of motion. In this case, the problem can be reduced to one 
DOF and is integrable ( Lidov & Ziglin 1976). This feature has been 
accounted for in many recent papers, moreover, the apparently sub- 
tle third-order perturbation to the Keplerian model, or the first order 
perturbation to the integrable quadrupole-order approximation may 
introduce qualitative changes of the dynamics. 

We calculated the secular expansion (Eq. [15]) up to the 10-th 
ordei0 One should be aware that by increasing the order of this 
expansion, we do not necessarily improve the approximation of the 
secular model of the real system, because this model is still lim- 
ited by the first order perturbation theory. In Section 3.2 we will 
examine the accuracy of the secular expansion in more details. 

Finally, the averaged relativistic correction possesses the same 
form as in coplanar case (Migaszewski & Gozdziewski 2008). 
Moreover, we include this perturbation only to the mutual inter- 
action of masses mo and m\ : 



3(A/t) 4 (Pt) 
c 2 L\G x 



+ const, 



(18) 



as it was explained above. 

Having the averaged model in hand, we may calculate the sec- 
ular frequencies of the inner companion, and compare them with 
the mean motion of the outer object (n 2 ). For the relativistic ad- 
vance of the inner periastron we have: 




2 This expansion is available on request in the form of a raw MATHE- 
MATICA input file; it will be also available on-line, after publishing this 
manuscript. 
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and for the apsidal motion forced by the mutual interaction of the 
inner and outer companion (in the quadrupole approximation): 



/l,ir 



n 2 




/Afi m2 a 3/2 



1 



=Ai- 



m\ 2 
-or 



1 



m " (l-ej) 2 



A 2 , 



where Ai^ are the following functions of the geometric elements: 

Ai = ^(1 - e?) [(3C 7 2 - 1) - 5(C 7 2 - l)cos2co! 

+ \c] [(2 + 3e 2 )-5e 2 cos2coi 
8 L 

A 2 = \ci [(2 + 3e 2 )-5e 2 cos2coi 
8 L 



Assuming now that fi\ 2 ~ & 2 mo (the central mass dominates), we 
may obtain the following estimates of the secular frequencies in 
terms of the characteristic units in our model: the relativistic fre- 
quency relative to n 2 is 



n 2 V lm 



0.2au\ /0.04au\ 



a\ 



3/2 



1 



a / 



while the mutually forced frequency relative to n 2 : 



fl,n 



ll 2 



;x 10" 



m 2 
10m, 



lm© 
m 



a \3/2 



(—) 
V 0.047 



(i-^) 3/ V 



=Ai 



+ 2 x 10- 



rn i 
100m, 



( lm & \ ( a 



J\m ) WoJ (I -el) 



;A 2 . 



These frequencies may be compared with the tidal apsidal fre- 
quency induced by the primary and the inner body bulge (Mi- 
gaszewski & Gozdziewski 2010, in preparation): 



/i,tid 
n 2 



4x 10" 



111 j 



+ 7x10 



lmQ 

lOOrnj J \ mo 
mo \ ( lOOrnj 



go 

1% 

5 



kL,0 

0.03 



lmQ ) \ m\ J \2Rj 

0.2m\ 5 ( 0.04 \ 3/2 l+3e 2 /2 + e 4 / 8 
J (I-* 2 ) 5 



0.15 



where Ro is the stellar radius, R\ is the radius of the inner sec- 
ondary and are tidal Love numbers of these bodies. Let 
us choose a reference model through setting characteristic parame- 
ters of a\ ~ 0.2 au, a ~ 0.04, mo ~ lmQ, mi ~ lOOrnj, Ro ~ IRq , 
/?i ~ 2/?j. Assuming that the bodies are modeled by polytropes with 
indices of 3 and 2, respectively, we compute their Love numbers, 
~ 0.03 for the primary, and ~ 0.15 for the inner secondary. Then 
setting e\ ~ 0, we obtain that the relativistic frequency is compara- 
ble with the mutual, Newtonian frequency, while the mean motion 
of the outer secondary is orders of magnitude larger than both of 
them (hence no mixed resonances are possible). Simultaneously, 
the assumptions of the averaging principle are well fulfilled. This 
guarantees that the evolution of the mean (secular) system closely 
follows the real configuration over the time scale of order ~ 1 /e, 
where e is the small parameter of the perturbation. 

Under the same assumptions, the tidal frequency is orders 
of magnitude smaller than the leading frequencies of the mutual 
(Newtonian) and relativistic corrections. This means that the tidal 
effect is negligible, indeed, as far as the model parameters do not 
strongly deviate from the characteristic values, as defined above. 



3.1 A global, 2-dim representation of the phase space 

Because the general planetary Af-body problem is very complex, we 
restrict the further analysis to its simplest, non-trivial case of three 
bodies ("non-trivial" in the sense of its non-integrability). We shall 
consider configurations of the host star and two planets or C-type 
systems comprising of a binary and a more distant body (a planet). 

We recall that the secular Hamiltonian i^ ec of the three body 
problem does not depend on M\,M 2 , therefore the conjugate ac- 
tions (L\,Li) are constant. The Hamiltonian 9l written in the 
Laplace reference frame depends on A£l = Cl 2 — Cl\ = ±7t only, 
not on Q.\ and ^2 separately. Hence, the following canonical trans- 
formation (e.g., Michtchenko et al. 2006 ): 



(©i,Gi) 
(co 2 ,G 2 ) 



(0)i, Gi) 
(co 2 ,G 2 ) 
(6i = £(Qi 

(e 2 = £(fli 



+ ^2),/l : 
-^2),^2 : 



--H1+H2) 
--H1-H2) 



(19) 



removes ^1,^2 from the secular Hamiltonian. After this transfor- 
mation it does not depend on 81, therefore J\ = |C| = C = const, 
where C is the total angular momentum of the system. Moreover, 
62 = ±7i/2 = const (after the Jacobi reduction of nodes) and J 2 
may be expressed as a function of G\ , G 2 and J\ in the following 
form: 

J 2 = (G 2 l -G 2 2 )/J l . 

Therefore, for constant values of the angular momentum J\=C and 
the secular energy i^ ec , the secular dynamics are reduced to two 
DOF Hamiltonian system. Instead of the total angular momentum 
C, we shall use the so called Angular Momentum Deficit {AMD) 
( |Laskar|2000| ): 

AMD = Li+L 2 -C, 

or its normalized value of A = AMD/{L\ +L 2 ), 9L G [0, 1] ( |Mi-| 
gaszewski & Gozdziewski 2009a). Because L\, L 2 and C are inte- 
grals of the secular system, the relativistic correction, Eq. 18 does 
not change J^, and the DOF number does not change. Because 
(JHqsl) depends on G\ only, thus it affects only the temporal evo- 
lution of coi. We note here that the perturbation induced by the 
quadrupole moment of the star, which was discussed ( Migaszewski 
& Gozdziewski 2009b ) in the coplanar case, also depends on H\ in 
the 3D problem, i.e. on the orbital inclination to the equatorial plane 
of the star. This perturbation introduces an additional frequency to 
£2i and then AO. is not constant anymore. It means that we could 
not perform the reduction of nodes and the secular problem would 
have three DOF. This also means that the Laplace reference frame, 
defined in terms of the total orbital angular momentum, does not 
possess any constant orientation in space. Being aware of this prob- 
lem, we do not consider the dynamical flattening of the star and/or 
of the innermost planet. The two DOF model is then less general but 
the dynamics are better tractable, thanks to the geometrical tools, 
which can be applied to study this basic, low-dimensional problem. 

If we fix the secular Hamiltonian in the form of ^ ec = 
^ ec (Gi,G2,C0i,C02), men ma y be considered as a free param- 
eter of the secular model. Moreover, the phase space is four- 
dimensional, and to represent the phase space of the system glob- 
ally in terms of two-dimensional sections, which are easy to visual- 
ize, we follow a concept of the representative plane of initial condi- 
tions ( Michtchenko & Malhotra 2004 ), the E-plane from hereafter. 
The L -plane may be chosen in different ways, although all repre- 
sentations may be fixed and defined through the following condi- 
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tions: 



:0, 



:0. 



(20) 



the mean Hamiltonians, derived through the analytic ("A") and nu- 
merical ("N") algorithms may be defined as follows: 



These conditions imply that all phase trajectories of the secular sys 



tern cross the Z-plane (Michtchenko et al. 2006; Libert & Henrard 



2007), see also our explanation in (Migaszewski & Gozdziewski 
2009a). In accord with the symmetries in the secular 3D model, the 
solutions to these equations are the following four pairs of angles: 

(CO! , co 2 ) = {(0, 0); (0, ±ti); (±w/2, ±n/2); (±7t/2, ^%/2)}> (21) 

that also define four distinct quarters of the Z-plane, numbered with 
Roman numbers II, I, IV, and III, respectively, see (Michtchenko 
et al. 2006) for details. Let us note that no other combinations of the 
angles are permitted by Eqs.[20| This feature of the secular system 
flows from the symmetry of the secular Hamiltonian with respect 
to the apsidal lines of the mean orbits, and may be also justified 
by the explicit form of the equations of motion derived from the 
expansion of the perturbing Hamiltonian, see ( Michtchenko et al.| 
2006, Migaszewski & Gozdziewski 2009a) for details. 

In this sense, the Z-plane may be thought as an analogue of 
the Poincare cross section. The conditions fixing the characteristic 
plane may be also rewritten as follows: 

coscoi = cos CO2 = U sincoi = sinc02 = 0. 

Further, we shall use three, basically equivalent geometric repre- 
sentations of the E-plane, which cover certain combinations of the 
quarters (the solution pairs of the pericenter arguments): 

• the ^-plane is defined with coscoi = cosC02 = 0, and 

fPs = {x = eismC0i,y = 6? 2 sinC02,ei 6 [0,1), e 2 € [0,1)}, (22) 

• the ^c-plane is defined with sincoi = sinC02 = 0, and 

fP c = {x = eicosC0i,y = £?2COSC02,ei e [0,1), e 2 G [0,1)}, (23) 

• and, finally, two incarnations of the 5-plane: 

S = {x = e\ cosAG5,y = £?2cos2c0i,ei e [0, \),e 2 £ [0, 1)}, (24) 

5' = {x = e\ cos2c0i,y = <?2COsAG5,£i G [0, 1),<?2 ^ [0, 1)}, (25) 



that was defined originally in (Michtchenko et al. 2006). It can be 
shown, that the and ^c-planes carry out the same information 
as the 5-plane. However, the 5-plane has a discontinuity along the 
x-axis, and the former representations are sometimes more conve- 
nient to the analysis of solutions of the secular system. 



3.2 A test of the analytic model 

We left a test of the accuracy of the secular expansion to this end, 
because the introduced E-planes are useful to illustrate the results 
of this test in a global manner. We select initial conditions in the 
S -plane, and the secular energy computed with the help of the an- 
alytic expansion is compared with the results of numerical averag- 
ing developed in (Gronchi & Milani 1998, Michtchenko & Mal- 
|hotra|2004| ), which are exact up to the numerical quadrature error. 
We consider the non-relativistic case only, because the secular rel- 
ativistic correction is exact (with the first non-zero post-Newtonian 
term included), and it does not influence the precision of analytic 
formulae. Figure [I] shows the levels of ^ ec , marked with solid 
curves in the 5-plane. Each panel is for a different order of the 
analytic approximation. The relative difference between values of 



■^sec 



(26) 



where / is the order of the analytic expansion in a. The results 
of this comparison are illustrated in Fig. [T] that shows the lev- 
els of A/ computed for a hierarchical system with a = 0.04 and 
ju = m\/ni2 = 5. The quadrupole-level model reproduces the sec- 
ular Hamiltonian as the even function with respect to both x = 
ei{sin,cos}cOi and y = £2{sin,cos}c02. The higher order approx- 
imations of 5^ ec broke this symmetry. We have shown in ( |Mi-| 
gaszewski & Gozdziewski 2009a) that the shape of 9-Q &c signifi- 
cantly depends on a. This is more important in the spatial prob- 
lem, because even for relatively small a, the quadrupole model dis- 
torts the structure of the phase-space (see Sect. 4.1). An inspection 
of Fig. [T] reveals that the octupole model reconstructs the secular 
Hamiltonian and its shape in the 5-plane very well, because the 
largest deviations A3 ~ 10 -3 appear only for e\ ~ 1. In other parts 
of the representative plane, A3 ~ 10 -5 . The high-order expansions 
are obviously even more precise. Following estimates of the secu- 
lar frequencies in the relevant parameter ranges (see Sect. 3), the 
averaging principle assures us that the orbital evolution of the sec- 
ular system follows closely the "real" orbits. Hence, we may quite 
safely skip a comparison of the results from both approaches by 
direct numerical integrations. 

The octupole model introduces the first order perturbation to 
the integrable quadrupole model, hence, because it is very precise 
in the range of small a, we further focus on this most simple, non- 
trivial case. 



3.3 Equilibria in the secular 3D problem 

In this work, we focus on the simplest class of solutions that are 
the equilibria (or stationary solutions). These solutions imply the 
basic structure of the phase space. By determining their stability, 
we may derive the local structure of neighboring phase trajecto- 
ries through relatively simple analysis. The equilibria of the non- 
resonant system in terms of the quadrupole approximations were 
studied in the past, (e.g., Krasinski i|1974[|Lidov & Ziglin 1976]). 
In our earlier work ( Migaszewski & Gozdziewski 2009a), we clas- 
sified families of equilibria emerging in the two-planet, non-planar 
problem, in terms of the quasi- analytic averaging, and basically ex- 
act secular model. We studied a few families of equilibria known 
in the literature, e.g., the zero-eccentricity solutions and the Lidov- 
Kozai resonance. We also found new families of these equilibria, 
in particular, the so called chained orbits solution that could be 
hardly derived with the perturbative approach, although the results 
in |Gronchi & M ilani (1998) might be applied here. This work fo- 
cus on the planetary regime of parameters ^,a, the mass ratio \i 
was restricted to the range of [0.1,2] and a e [0.1,0.667]. More- 
over, we explored the whole permitted range of G [0,1]. Yet we 
learned that the semi-analytic approach has serious disadvantages. 
All calculations, including the integrations of the equations of mo- 
tions, and the stability analysis must be performed with the help of 
numerical algorithms. This may introduce large errors and hinders 
the analytic, qualitative analysis of the problem. 

In this section, we extend the study of stationary solutions in 
( Migaszewski & Gozdziewski 2009a) to a wider range of the mass 
ratio, ju > 2. Simultaneously, we consider smaller a < 0.1. That 
makes it possible to apply the analytic model described and tested 
in Sect. 3. Because planets emerge most likely from remnants of a 
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Figure 2. Levels of ^ ec at the ^-plane calculated for a = 0.01 and \jl = 20. Panels from a) to f) are for different values of A = 
{0.08,0.12,0.16,0.19,0.2,0.25}, respectively. These values of A imply the mutual inclination at the origin i = {49°, 61°, 70°, 77°, 79°, 89°}, respec- 
tively. The intensity of shaded ares encodes the mutual inclination in prescribed ranges, darker shade is for larger i mut . The inclinations ranges are 
z'mut = (35° , 45° ) , (40° , 55° ) , (50° , 65° ) , (60° , 75° ) , (60° , 75° ) , (70° , 85° ) in subsequent panels. 



thin protoplanetary disk, we also restrict our attention to mutual in- 
clinations up to z' mut ~ 90° (direct orbits). In that range, we may find 
families of equilibria related to the zero-eccentricity orbits, and the 
LK resonance classified in ( Migaszewski & Gozdziewski 2009a) 
as solutions of family IVa, accompanied by families Ilia, Illb, and 
IVb+. Our analysis is also restricted to the initial conditions in quar- 
ters III and IV of the S -plane, in which these solutions may only 
"reside". This region of the parameter space has been studied (e.g., 
Krasinskii |1974| ) with the quadrupole-level model which is some- 
how the next to trivial, non-interacting Keplerian approximation of 
the three-body orbits. However, as we show below, this approxima- 
tion may introduce artifacts due to generally non-realistic symme- 
try of the secular Hamiltonian. Obviously, to avoid this problem, 



higher order expansions are required. We also attempt to extend the 
results of |Ford et al. 1(2000) , who applied the octupole-level theory 
to the analysis of hierarchical planetary systems with very small 
a -0.01. 

To show the generic properties of the relevant families of equi- 
libria, we begin with an example that is illustrated in Figure[2] Pan- 
els in this figure show levels of ^ ec at the fP^-plane for a = 0.01 
and fi = 20 and a few different values of Jl. Panel [2^ was derived 
for 2L = 0.08, and it reveals the zero-eccentricity equilibrium. For 
larger A = 0.12 (Fig.j^) the fi sure-eisht structure of Lidov-Kozai 
equilibrium appears, and is labeled with IVa. We recall, that the 
mutual inclination of circular orbits that corresponds to the LK 
bifurcation will be called the critical inclination, z' crit although, in 
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-0.5 -0.25 0.25 e^ino^ 0.5 



Figure 3. Families of stationary solutions (IVa, Ilia, Illb, IVb+). Semi- 
major axes ratio a = 0.01, mass ratio \x = 20. Dark dots are for stable equi- 
libria, grey dots are for unstable equilibria. Crossed and dotted circles mark 
bifurcations. Small arrows shows the direction of particular stationary solu- 
tion with increasing SL. Green dots are for the positions of equilibria for par- 
ticular values of J^ = 0. 12,0. 16,0. 19, 0.2, 0.25 (the energy levels are shown 
in Fig.|2j. Solutions are for classic model, and were obtained with the help 
of the octupole theory. 



general, any such value of the mutual inclination that leads to a 
bifurcation of equilibria has the sense of being "critical" ( |Krasin-| 
|skii|1972|[l974] ). For larger value of 2L = 0.16 (Fig. [2}:) the LKR 
"moves" towards larger e\ (simultaneously, £2 ~ ) and a new sad- 
dle point appears. We call this solution a member of family Ilia. 
For larger Si = 0.19 (Fig. [2jl), these structures still expand, and 
for Si = 0.2 (Fig. [2^) two new equilibria emerge: one is associated 
with quasi-elliptic point (family Illb from hereafter), and a saddle 
of family IVb-h Apparently, it emerges from a point near the IVa 
solution in the fP^-plane, nevertheless it does not correspond to a 
bifurcation of this equilibrium, see (Migaszewski & Gozdziewski 
2009a). Solution IVa moves towards larger e^, see Fig.|2p. 

The parametric paths of these equilibria in terms of Sk, may 
be depicted in the ^-plane (Fig. [3]). Black and grey curves are for 
the stable and unstable equilibria, respectively. The relevant fami- 
lies of equilibria are labeled with Roman numbers and Latin letters. 
The direction of "motion" of particular solutions along the J^-axis 
is marked with arrows. For a reference, positions of the equilibria 
for a few discrete values of SA. = 0.12,0.16,0.19,0.2,0.25 (corre- 
sponding to subsequent panels in Fig. [2]), are marked with green 
dots and labeled. Following a particular evolution path of equilib- 
rium IVa, we see that it appears for Si. ~ 0. 1 through a bifurcation 
of the origin. When SA increases, this solution moves along £2 ^ 
towards larger e\. For SA ~ 0.17, it reaches the maximal e\ ~ 0.4, 
and it turns back, towards smaller e\ with simultaneous increase of 
^2- When Si. increases even more, this solution reaches the border 
of convergence of the analytic expansion. We call that border the 
anti- collision line [see (Migaszewski & Gozdziewski 2008)]. 

The parametric evolution of equilibria in quadrant III is more 
complex. The first non- trivial stationary solution appears for SA ~ 
0.15, which is unstable, saddle point Ilia. It evolves along e\ ~ 
and then £2 increases to large values. For SA ~ 0.195, two new solu- 
tions Illb and IVb+ emerge from the elliptic structure related to the 
LKR (see Fig. [2^). One of them is stable (solution IVb+) and the 
another one is unstable (solution Illb). They appear around (e\ ~ 
0.45,^2 ~ 0.15) . When increases, the solution IVb+ moves 
towards e\ —> 1 and e^ — > 0. Simultaneously, equilibrium Illb is 
shifted towards increasing £2 and decreasing e\. Then also Ilia in- 



creases £2 and leaves off the e\=Q axis. Finally, for Si ~ 0.26, 
equilibria Ilia and Illb merge at one point {e\ ~ 0.13,^2 ~ 0-8). 



4 THE SECULAR CHAOS 

Following (Michtchenko et al. 2006 ), we shown in (Migaszewski 
& Gozdziewski 2009a} that the averaged 3D model may involve 
strongly chaotic motions, in the secular time-scale. To study in de- 
tails, how the model parameters influence the structure of the S- 
plane, and how it relates to the long-term chaotic phenomena, we 
apply the fast indicator approach. Among many numerical tools 
of this kind, we choose the so called coefficient of the diffusion of 
fundamental frequencies a introduced by (Laskar 1990 ). To check 
whether a phase space trajectory of a quasi-integrable Hamiltonian 
system is regular (quasi-periodic) or irregular (chaotic), one inte- 
grates the equations of motion over two subsequent intervals of 
time, e.g., [0,T] and [T,2T]. Next, we resolve the frequencies in 
the discrete orbital signal with the help of refined FFT analysis 
(Laskar 1990 ), obtaining two estimates of a given frequency, say 
fr and fxr, over these two intervals of time. The coefficient of the 
diffusion of the fundamental frequency is then defined through: 



Clearly, if the signal does not change over time, a ~ and this 
means that the phase trajectory is quasi-periodic (stable). If a is 
significantly different from 0, the trajectory is chaotic and regarded 
unstable. In our calcul ations, we used a variant of the frequency 
analysis developed by |Sidlichovsky & Nesvorny|l996| , which is 
called the Frequency Modified Fourier Transform (FMFT). We also 
used publicly available code of the FMFT algorithm, kindly pro- 
vided by David Nesvorny on his personal web-pag^] 

Because the secular evolution is associated with 00/ angles, 
we compute the a coefficient on the basis of complex time- series 
{Gift) expiGfy(f)}, . . . i = 1,2, where i is the imaginary unit. In this 
signal, the osculating eccentricity and pericenter argument for each 
planet represent rescaled canonical action-angle variables. Hence, 
resolving its Fourier components, we may determine the leading 
amplitudes (proper eccentricities) and the fundamental frequencies 
of pericenter angles. 

Next, we did massive integrations of the secular equations of 
motion. The initial conditions were selected at the grid of 200 x 100 
data points of the ^-plane. At each point of the dynamical map, we 
integrated the secular trajectory over ~ 10 4 secular periods with 
respect to the smaller frequency (typically, one of the fundamental 
frequencies is much larger than the other one). Having the com- 
puted phase trajectory, we then find an estimate a, as well as the 
maximal eccentricities attained by both orbits during the integra- 
tion time (the so called max^ indicator) as well as the amplitude 
of variation of the mutual inclination A/ = (max/ mut — min/ mut ) at- 
tained during the integration time. These geometrical characteris- 
tics are very useful to understand the source of instabilities indi- 
cated and detected by the diffusion coefficient a. 

Figures |4}j6] illustrate the results derived for the same system 
that energy levels are shown in Fig. [2] Subsequent figures are for 
SA = 0.12,0.16,0.2, respectively (see Figs. [2]3,c,e for the respec- 
tive levels of ^4c)- (We note, that for SA = 0.08 the view of the 
phase space is basically very simple and the motions are regular 

3 http://www.boulder.swri.edu/~davidn/ 
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Figure 4. Dynamical maps for the classic (point-mass Newtonian) model, shown in the IPs-plane. Panels from a) to d) are for the coefficient of the diffusion 
of fundamental frequencies (a), maximal eccentricity of the inner and outer planet (maxei ,max^2), respectively, and the amplitude of variation of the mutual 
inclination (Ai). Dots mark areas of librations of C0i (panel b), CO2 (panel c), and AG5 (panel d). Stationary solutions are marked with circles: dotted circles are 
for stable equilibria (corresponding to the maximum of the secular Hamiltonian), crossed circles are for unstable equilibria, empty circles are for the linearly 
stable equilibria. These solutions are labeled, in accord with ( Migaszewski & Gozdziewski 2009a). Parameters of the system: a = 0.01,/u = 20,J% = 0.12, see 
also Fig.[2j). 



0.2 



0.4 



max ei 



Illal 



V 



0.9 

sT 



0.3 



IVa 




0.9 



0.3 



-0.8 



-0.4 



10 



0.4 



Ai [deg] 20 



d) 




A=0.16 



0.9 



- 0.3 



0.4 e^incOj 



-0.8 



-0.4 



0.4 e^incOj 



Figure 5. Dynamical maps for the Newtonian, point mass model in the ^s-plane, = 0.16. See the caption to Fig.[4]for more details, and also Fig.[2j:. 



everywhere). The right-hand panel of each figure is for a. The dy- 
namical character of the phase- space trajectories is color-coded: 
black color means quasi-regular evolution of the secular system 
(a ~ 10 -6 -10 -8 ), and yellow colour is for strongly chaotic mo- 
tions (a ^ 10 -4 ). The a-maps reveals that almost the whole phase 
space is filled up with regular motions, and chaos appears only in 
some small regions in the ^-plane. For a reference, the coordinates 



of equilibria IVa, Ilia, Illb, IVb+ and (the equilibrium at the ori- 
gin) are marked with circles. Dotted circle is for Lyapunov stable 
solution IVa, crossed circles are for unstable solutions 0, Ilia and 
Illb, respectively. Equilibrium IVb-i- is linearly stable and is marked 
with open circle. Clearly, trajectories close to unstable equilibria 
Ilia and are chaotic. Let us note that unstable equilibrium Ilia lies 
in a very narrow, chaotic zone too. 
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Figure 6. Dynamical maps for the Newtonian, point-mass model, in the fP^-plane, Si = 0.20. See the caption to Fig.[4]for more details, also Fig. [2^. 



ble, similarly to the LKR in the restricted problem. This dynamical 
phenomenon may explain a large eccentricity of the binary, when 
a distant, low-mass and dark companion cannot be detected due to 
very long orbital period. 

The dynamical maps in Figs. |4jj6]re veal some zones, in which 
the outer eccentricity is strongly amplified. This eccentricity is ob- 
viously constant in terms of the quadrupole model. This feature of 
the octupole model shows that the small perturbation may cause 
extended, geometric changes of the mean orbits. The amplifica- 
tion of e2, with simultaneously almost constant relative inclination, 
Ai ~ const (see Figs. |6j;,d) seems appear due to the bifurcation of 
equilibria IVb+ and Illb for SA ~ 0. 195. 



4.1 A model explaining the secular chaos 

The origin of chaotic secular dynamics may be explained by the 
presence of separatrices, which encompass different types of mo- 
tions, librations and rotations of angles ©i, CO2 and AG5. A classifi- 
cation of these libration modes in the secular problem of two plan- 
ets modes has been introduced in ( Michtchenko et al. 2006 ). The 
separatrices appear due to equilibria and periodic solutions in the 
integrable, or close to integrable secular models, which might be 
understood as the quadrupole approximation and/or the co-planar 
configuration. 

Here, we found a simple explanation of the mechanism gen- 
erating chaos, which is in fact the same as in the perturbed pen- 
dulum. To show this, let us recall that the expansion of SH^ QC to 
the second order in a is the celebrated integrable quadrupole-level 
model. The energy levels of this model are illustrated in Fig. [7^. 
Because e^ is the integral of motion, the representative plane may 
be constructed similarly to the integrable co-planar problem, S /r = 
{e\ cos2c0i x £2}- In the region marked in green colour, any con- 
stant level of e2 crosses the energy curve in two turning points limit- 
ing the range of variation of e\ . If the e^ level is tangent to the given 
level of the energy, the dynamics must be then confined to fixed e\, 
hence we obtain an equilibrium point (stable or unstable). A set of 



Besides dynamical maps for a, Figs. |4jj6] illustrate geometric 
characteristics or orbits in the ^-plane, in terms of max e\£ indica- 
tor for the inner and the outer orbit, respectively, and for the ampli- 
tude of variation of the mutual inclination, Ai. In all these dynam- 
ical maps, the green dashed contours mark the mutual inclination 
equal to the critical inclination of the LKR bifurcation ( Krasinskii 
[T972l[T974| ). Also regions, in which the secular angles C0i, CO2, and 
AU5 oscillate, are shown: the gray/white/black dots surrounded by 
appropriate curves bound initial conditions that lead to librations of 
COi around ±7t/2 in the max e\ -maps, librations of 002 around ±7i/2 
in the max £2 -maps, and librations AG5 around 0,71 in the A/-maps. 
We note that 00 1 librates in almost all regular trajectories with the 
initial / mut > / C rit • If this condition of the Lidov-Kozai mechanism 
is fulfilled then COi librates around ±7c/2. In chaotic zones of the 
5-plane, these angles do not oscillate, even if condition / mut > z' crit 
holds true. Angle 002 librates around 7i/2 only when e^ ~ 0, for 
SA = 0. 12, 0. 16. This region extends significantly for larger Si = 0.2 
. In such a case, there are three zones in which both angles COi and 
CO2 librate: a region connected to equilibrium IVa, the neighborhood 
of linearly stable solution IVb+, and a region of small e\ (though 
there is no stationary solution associated with these librations). The 
later area is surrounded by strongly chaotic motions shown in the 
G-plane. 

Clearly, the equilibria permitted by relatively large SA affect 
the secular dynamics, which become very complex. This may be 
better seen in the max^-maps. Subsequent panels in Figs.[4}j6]show 
that the inner eccentricity may be strongly excited if z mut > z crit . For 
relatively small $1 = 0.12, the inner orbit, which is initially quasi- 
circular, becomes moderately eccentric, with max^i ~ 0.3. Note 
that in this specific case, the inner body is 20 times more mas- 
sive then the outer companion. Simultaneously, as Fig.|4]shows, the 
outer eccentricity e^ is not amplified. We may conclude that in the 
regime of the outer LKR, when the inner body is much more mas- 
sive than the outer body, the secular evolution may lead to strong 
perturbations of the inner orbit. Hence, even the mass hierarchy is 
reversed, a strong excitation of the inner eccentricity is still possi- 
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these equilibria for increasing, fixed values of £2 is marked by two 
thick curves that meet around of {e\ = 0,^2 ~ 0.4). This is the bi- 
furcation point, at which two families of equilibria emerge — the 
stable branch of the LKR and the unstable origin (e\ = 0, = 0). 

For a given value of constant e2, we may then find the max- 
imum range of e\, which corresponds to librations of the canoni- 
cal angle coi , and this value determines the position of separatrix, 
see Fig. for the separatrix shown in the phase diagram in the 
{e\ cos2coi x e\ sin 2coi} -plane, corresponding to the energy level 
marked by the blue curve on the levels diagram (Fig.[7^). Collecting 
such points along increasing e2, we construct the red dashed curve 
in Fig. [7^, showing the separatrix region globally in the whole 
-range (see the shaded area in Fig. [7^). The dashed line that marks 
the unstable origin, might be also understood as the other branch of 
the global separatrix, which, for any value of e^, corresponds to the 
equilibrium point. 

Now, considering the perturbed quadrupole-model in terms of 
the octupole-level secular expansion, we may expect that chaotic 
motions will appear close to the separatrix, due to the perturbation. 
Indeed, this is illustrated in two panels of Fig. [8] The left panel is 
for the frequency diffusion a computed for initial conditions in the 
S /r = {e\ cos2coi x ^j-plane. This plane shows the energy levels 
of the octupole model, with stable (thick, solid curve) and unstable 
(thick, dashed line) equilibria in the quadrupole model over-plotted. 
The red curve is for the global separatrix of the LK resonance in the 
quadrupole model, and the green region is for the librations of COi 
in the integrable case. Clearly, chaotic orbits are found along the 
branch of unstable equilibria, as well as close to the red separatrix 
curve. This may be better seen in Poincare cross section computed 
for a fixed energy level marked with blue curve that is shown in 
the right-hand panel of Fig. [SJd. This panel in fact comprises two 
sections AG5 = n in the {e\ cos2coi x e\ sin2coi }-plane, for AG5 > 
(the left half-plane), and for AG5 < (the right half-plane). In place 
of the separatrix curves, a wide chaotic regions appear. 

A careful inspection of Fig. [8^ reveals additional narrow 
banana- shaped chaotic structures. They are related to the separatrix 
of the librations of angle AG5 around n and simultaneous circulation 
of angle COi, which is classified as mode 2 in ( Michtchenko et al. 
2006 see their Fig. [3]) . Further, the plots of fundamental frequen- 
cies computed along the x-axis of the S '-plane in Fig. [8^, which are 
shown in Fig. [9] reveal that indeed, in this region the fundamental 
frequencies decrease to 0, indicating the separatrix region. 

This example shows clearly that apparently very small per- 
turbation (recall that a = 0.01) to the integrable model introduces 
extended chaotic behavior and qualitative change of the secular dy- 
namics of the system. 

This interpretation is helpful to understand the source of 
chaotic motions shown in all dynamical maps (e.g., Figs. |4|{6} - 
usually, they appear close to the separatrices associated with unsta- 
ble equilibria or resonances of the secular angles (unstable periodic 
orbits). 

Figure [l0| shows the temporal evolution of the mean orbits 
with initial conditions close to the origin in Fig. [4] (see its caption 
for the details). The eccentricity of the inner orbit is significantly 
perturbed while the outer orbit is almost unchanged. The mutual 
inclination evolves in anti-phase with e\, and COi librates around 
7i/2. These are typical features of the Lidov-Kozai resonance, still 
we recall that in this instance, the smaller outer body forces the LK 
cycles on the orbits of much more (20 times) massive than the inner 
component of the binary. 
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Figure 10. The mean orbits with the following initial condition: mo = 
1 m , m\ = 200 mj, mi = 10 mj, a\ = 0.3 au, ai = 30 au, e\ = 0.001, 
e 2 = 0.001, COi = 7T/2, co 2 = 7T/2, z mut = 59°.64, A = 0.12. Panels from the 
top to the bottom present the time-evolution of <?i,<?2,*mut, Wi,C02, respec- 
tively. 



5 EQUILIBRIA IN THE CLASSIC OCTUPOLE MODEL 

The dynamical maps and their analysis show that the equilibria 
constrain the secular dynamics of the perturbed model. In terms 
of the Newtonian, point-mass formulation, the stationary solutions 
depend on parameters a and ju. Limiting our survey to / mut < n/2 
(direct orbits), we perform a parametric survey of the equilibria 
in terms of the octupole expansion, in such a manner that a com- 
parison with the results derived for the more general, relativistic 
model will be possible. Here, we consider a more extended range 
of mass ratio ju, covering a transition from the planetary regime 
(small jjl) to the circumbinary case (large ju ) than in (Migaszewski 

6 Gozdziewski 2009a). Yet the assumption of small a makes it 
possible to use the analytic formulation of the secular Hamiltonian, 
which is very precise in terms of the octupole-level approximation, 
instead of the numerical approach in (Migaszewski & Gozdziewski 
|2009aT >. 

The parameter dependence of the equilibria is illustrated in 
Fig.[TT]that shows the fP^ -plane (note that due to the symmetry, only 
the upper half-plane is illustrated, see also Fig. [3}. We set a = 0.04, 
jd e [1.5,3,5, 10,20,25,50], and then positions of the equilibria may 
be traced along increasing J3, which might be understood as the 
natural curve parameter. 

In fact, instead of \x, we choose a new parameter (3 (see |KrasirT| 
|skii|1974[|Migaszewski & Gozdziewski|2010| > 

p(A/,a)=Li/L 2 -^\/oc, 

that better reflects the dependence of the dynamics on both \i and a 
than on one of these parameters itself (we did a few numerical tests 
for different a which confirm this scaling very well). Actions L; 
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Figure 7. The left panel is for the levels of the secular energy (thin curves) depicted in the representative plane of the quadrupole problem, 5" = {e\ cos2coi x 
<?2 }-plane, a = 0.01 , /a = 20, A = 0. 12. Recall that this problem is integrable and <?2 = constant is a parameter. The black, thick curve marks the LKR equilibrium 
for different energies and ei integral. The red, dotted curve marks the separatrix between librations and rotations of C0i for fixed values of integrals. The shaded 
(green) zone indicates librations of C0i around 7i/2. The right hand panel is for the phase diagram in the {e\ cos2c0i x e\ sin 2c0i} -plane and a fixed energy 
level marked by blue curve in panel a. The separatrix and the region of librations of C0i are also marked with green colour. 
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Figure 8. The left panel is for the secular energy levels (thin curves) shown in the S'-plane of the octupole model, a = 0.01, \x = 20, A = 0. 12. Shaded areas 
are for initial conditions leading to chaotic evolution of the secular model. The red curve is for the separatrix of the LKR resonance in the quadrupole model, 
the black thick curve is for the equilibria in this model. The right hand panel is for the Poincare cross-section AG3 = n in the {e\ cos2o)i x e\ sin2o)i }-plane. 
It corresponds to a fixed energy level marked with the blue curve in panel a. The separatrix of the quadrupole problem and a region of librating C0i are also 
marked. 
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Figure 9. The spectrum of fundamental frequencies computed along the x-axis of Fig. [8] for the same set of parameters. Abscissas, in which the frequencies 
approach zero, indicate the separatrices of secular resonances of the respective angles, C0i (panel a), and AG3 (panel b). Parameters are: mo = 1 awq, m\ = 40mj, 
ni2 = 2mj, a\ =0.1 au, a^ = 10 au. 
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reflect the angular momentum partition between both components. 
If (3 rsj 1, (p. rsj 5), both secondaries are dynamically equivalent, if 
(3 > 1 {]d > 5) then the inner body "dominates" dynamically in the 
system, and if (3 < 1 then the hierarchy is reversed. To illustrate the 
stability of the equilibria, the curves are marked with black dots are 
for Lyapunov stable solutions, and grey dots mean unstable solu- 
tions. The curves are labeled with both ju and (3 parameters. 

In the right-hand half-plane (quarter IV, coi = CO2 = 7t/2), for 
z'mut < ft/2, only one solution appears that is the LK resonance. For 
small \a, it evolves with increasing $1 along the axis of e 2 ~ be- 
tween e\ = (when the first LK bifurcation appears, see Fig.pjD) 
and e\ ~ 1. After the second bifurcation of the LKR (Fig. [2p), 
two solutions emerge. One of them is the LKR for z mut > n/2 (see 
Migaszewski & Gozdziewski 2009a this solution is not discussed 
here). The second new solution moves along the axis e\ ~ 1 up to 
large e 2 . 

For ju^2, the LKR does not reach e\ = \ but it "turns back", 
with decreasing e\ and increasing e 2 . For larger mass ratio, the 
maximal e\ is smaller. The families of LKR for ju = 1.5,3 cor- 
respond to systems with (3 < 1, hence the dynamical hierarchy is 
reversed, and the eccentricity of the inner, more massive body is 
forced by the outer companion. The position of the LKR family 
moves to the region of smaller e\ , and for large enough ju, this solu- 
tion is confined to e\ ~ 0-axis and tends to large e 2 . The direction of 
parametric evolution of this equilibrium, corresponding to increas- 
ing is marked with an arrow. We see that for all ju this direction 
is the same. 

The view of the left-hand half-plane is more complex. As we 
shown in (Migaszewski & Gozdziewski 2009a), the first solution 
appearing in quarter III of the ^-plane is unstable equilibrium Ilia 
emerging due to the second bifurcation of the origin (e\ = e 2 = 0) 
(marked with 1^ in our paper). Then, with increasing 5i , this so- 
lution moves along e\ ~ towards large e 2 . For some value of 
(e.g., - 0.195 for a = 0.01, /u = 20; see Fig.[3j solutions Illb and 
IVb+ appear in the same point of the phase space (this is illustrated 
with crossed circle in the left-hand half-plane of Fig. [3j. Solution 
IVb+ then moves along e\ — > 1 and e 2 — > 0. Simultaneously, solu- 
tion Illb evolves towards a point marked with dotted circle. Solu- 
tion Ilia also moves to that point and for some critical both solu- 
tions merge and vanish. Figure [TT]re veals that the parametric paths 
of equilibria depend on parameter \i (cc = 0.04 was fixed in this 
test). The path of the LKR becomes closer to e\ ~ 0-axis for larger 
mass ratio. Also families of stationary solutions that are present in 
quarter III, become closer to the origin. The range of eccentrici- 
ties corresponding to these equilibria is very small for /i ^ 25. We 
may also notice that the bifurcation points (crossed and dotted cir- 
cles, respectively) tend to each other with increasing mass ratio, 
which is consistent with the transition between the planetary and 
the circumbinary regime. For ju ^ 25, these two points are merged. 
In this particular case, the structure of equilibria in quarter III of 
the fP^-plane is more simple than in the general case because only 
one solution persists, a stable equilibrium corresponding to nearly 
circular outer orbit. 

For different a < 0.1, the general, global view of the fami- 
lies of equilibria is very similar to the results presented in Fig. [TT] 
In fact, as we noticed previously, the parametric evolution of the 
equilibria is reflected by parameter (3(^/,a), and basically does not 
depend on the individual values of a and \i. 



6 THE PN CORRECTION TO THE SECULAR MODEL 

The results presented in the previous section illustrate the already 
well known feature of the classic model (see, e.g., Michtchenko 
|& Malhotral|20041 |Michtchenko et al.| [2006). In the approxima- 
tion of small masses, the secular dynamics of the Newtonian 2- 
planet, hierarchical model depend on the semi-major axes ratio 
and planetary mass ratio, and not on individual system parame- 
ters, (a\,a 2 ,mi,m 2 ). In our works (Migaszewski & Gozdziewski 
2009b 2010), we shown that this feature is not preserved in the 
more general model, including relativistic, rotational and tidal cor- 
rections to the Newtonian point-mass interactions. In these settings, 
the secular dynamics depend on the individual semi-major axes, as 
well as individual planetary masses. Because the overall structure 
of the phase space is characterized by the equilibria, we now at- 
tempt to show that deviations between these equilibria in the classic 
and relativistic models become more important when the system di- 
mension scales down (a\,a 2 decrease when a is const), and masses 
m\,m 2 are smaller, when their ratio ju is kept constant. 

The differences between the two coplanar models manifest 
itself through the shapes and localization of stationary solutions 
and depend on the individual masses and semi-major axes \Mi-\ 
gaszewski & Gozdziewski 2009b ). Now we can observe the same 
feature in the spatial planetary system, corrected for the relativistic 
perturbation. Figure [12] presents families of equilibria in the same 
manner as Fig. [TT] (due to the symmetry, only the v-positive half- 
plane of &s -plane is presented). Families of stationary solutions in 
the classic model are drawn with blue and violet curves for sta- 
ble and unstable solutions, respectively. Solutions in the relativistic 
model are plotted with black (stable equilibria) and grey (unsta- 
ble equilibria) curves. In this experiment, we varied the individual 
masses of the secondary bodies, still keeping their ratio \i = 10 as a 
constant. The masses are changed between (m\ = 100 m h m 2 = 
10 mj) to (mi = 3 m f , m 2 = 0.3 im), and the primary mass is 
m = 1 m ■ • 

The results illustrated in Fig. [12] confirm that even for large 
secondary masses, the parametric curves of stationary solutions in 
the realm of the classic model depends weakly on the individual 
masses. Yet similarly to the co-planar case, curves of equilibria 
calculated with the relativistic corrections differ significantly from 
the results obtained for the classic model. The deviations between 
both models are more significant for smaller masses of the sec- 
ondary bodies. For instance, the family IVa moves to the range of 
much smaller e\. Solutions of this type exist even for very small 
m\ and m 2 , which is just not possible in the Newtonian model ( |Mi-| 
gaszewski & Gozdziewski 2009a ). When the mutual perturbations 
between secondaries are small enough, the critical inclination in the 
relativistic model, that leads to the LK bifurcation, becomes larger 
than 7i/2 . Thus / crit increases with decreasing masses. 

In quadrant III of the fP^-plane (002 = — ©1 = =b7c/2), the struc- 
ture of equilibria is even more complex. For some critical values of 
masses (m 2 ~ 2.07 im) the parametric paths divide into two parts. 
The top part, characterized by larger e 2 , comprises of two unstable 
equilibria emerging from one bifurcation point, which then meet in 
another bifurcation point. The bottom part of the equilibria curve 
represents a saddle point, that changes its stability, from unstable 
Ilia-type solution to the stable solution IVb+. This branch is very 
similar to equilibria in the classic system with large ju (or large 
(3 > 5) observed already in Fig. [TT] however this takes place for 
quite a different value of (3 = 2. 

The stationary solutions are determined by the shape of the 
secular Hamiltonian. Its levels are plotted in the ^-plane, Fig.[T3] 
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Figure 11. Families of stationary solutions in the ^-plane. The semi-major axes ratio a = 0.04, /a = { 1 .5, 3, 5, 10, 20, 25, 50} (each curve is labeled accordingly 
with the value of /u). Dark dots are for stable equilibria, grey dots are for unstable equilibria. Crossed and dotted circles mark the positions on the equilibria 
curves where solutions bifurcate (see the text for details). Small arrows show the directions of the evolution of particular equilibrium with increasing A. 
Parameter (3 = L\/Li (see the text). Stationary curves corresponding to (3 = 1 are drown with thicker lines. Stationary solutions were calculated for the 
octupole classic model. See also Fig. [3] 




Figure 12. Families of stationary solutions presented in the IPs-plane, calculated for a = 0.04, ,u = 10, ai = 0.2 au,«2 =5.0 au,rao = 1 m and varied 
m2 = 10,4,2.7,2.07, 1.85, 1.3, 1.0,0.3 mj. Equilibria in the classic model are compared with the equilibria in the relativistic model. The mass of the outer 
body m2 labels each particular curve. Black dots are for stable equilibria, grey dots are for unstable equilibria of the relativistic model. Equilibria of the classic 
model are plotted with blue and violet dots for stable and unstable solutions, respectively. Positions of equilibria were calculated with the help of the octupole 
theory. 



Each panel in this figure is calculated for the same parameters jli, 
cc, 2L (but in this case, particular values of masses and semi-major 
axes are varied). Let us recall, that Fig. [13^ shows the phase space 
calculated in classic model, while subsequent panels [T3jp,c,d are 
for relativistic model, and different masses and semi-major axes, 
labeled in subsequent panels. 

If the masses are large (Fig. pjjp), we can see four elliptic 



points separated by four saddles (let us recall that the -plane is 
redundant, and the energy levels are reflected with respect to the 
origin, thus in fact we have only two unique elliptic points and only 
two saddles). The elliptic points may be identified with solutions 
IVa (coi = co 2 = ±n/2) and Illb (a>i = -co 2 = ±rc/2). The saddles 
correspond to solutions Ilia (e\ ~ 0) and IVb+ ~ 0)- When the 
masses decrease (still, ju is kept constant), the structure surround- 
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ing solution Illb becomes smaller and moves towards solution Ilia. 
Simultaneously, the saddle point IVb+ tends to the origin. For the 
masses small enough, (Fig. [T2j, solutions Ilia and Illb merge and 
vanish, while IVb+ "falls" into the origin. 

Figures [T4}fT6| shows the dynamical maps for the relativis- 
tic model, and are constructed in the same manner as maps in 
Figs. [4}|6] Subsequent figures correspond to masses and & used 
to plot the energy levels in Figs. fT3fr ,c,d respectively: Fig. [14] is 
for classic model, Figs. [15] and |16| are for the relativistic model 
with rri2 = 10 im and mi = 1.85 mj, respectively. The mass ra- 
tio is n = 10 in both instances. The order of panels and symbol- 
coding of equilibria, as well as coding libration zones of the angles 
coi, 002, and AG5 are the same as in Figs.[4}j6] in particular, dot- 
ted, crossed and empty circles mark the Lyapunov stable, unstable 
and linearly stable equilibria, respectively. Clearly, the overall view 
of the phase space is different in all cases. The regions of chaotic 
motions (yellow areas in the a-maps) obtained for the classic and 
relativistic model are significantly different. Also the dependence 
on the masses of secondaries in the relativistic models is evident. 
The structure of chaotic/regular secular evolution is reflected in 
max -maps and through librational regions of coi and CO2. We 
note, that in the case of the regular solutions, coi librates around 
±n/2 when / mut > / crit . In some part of this region also CO2 librates 
around ±7t/2. 

Figure [TTJillustrates the temporal evolution for an initial con- 
dition written in the caption. The parameters of this model are 
the same as in Fig. [14] We chose the same initial eccentricities 
e\ = 0.45 and = 0.001, and integrate the secular equations of 
motion of classic (grey curves) and relativistic (black curves) mod- 
els. We note qualitative differences between both configurations. 
The classic model leads to much larger variations of the elements. 
Particularly, the outer eccentricity is strongly amplified, com- 
pared to the variations in the relativistic model. Still, this is not a 
rule. Inspecting the bottom row of Fig. ^] we can find regions in 
the fP^-plane, in which, for the same initial condition, max e^ be- 
comes larger in the relativistic model than in the Newtonian model. 
Also the secularly chaotic configurations appear in quite different 
zones of the phase space in both models. Moreover, the relativis- 
tic corrections may transform the regular evolution in the classic 
model into the chaotic evolution in the relativistic systems, and vice 
versa. Remarkably, the configuration illustrated in Fig.[l4]has very 
large masses m\ = 100 m, and mi = 10 m , . 



7 CONCLUSIONS 

In this work, we attempt to show that the global features of the sec- 
ular dynamics of the 3-D, non-resonant planetary system depend 
qualitatively on the apparently subtle relativistic corrections to the 
Newtonian gravity. A lesson, which we learned studying the co- 
planar case (Migaszewski & Gozdziewski 2009b) is that the non- 
Newtonian point-mass interactions might be very important for the 
global dynamics, because, in contrary to intuition one may have, 
the corrections to the Newtonian interactions might be not small, 
as compared to the mutual point-to-point gravity. The numerical 
analysis of this multi-parameter problem is complex, hence the un- 
derlining idea of this paper lies in a construction of possibly pre- 
cise analytical model. Essentially simple averaging of the perturb- 
ing Hamiltonian in (Migaszewski & Gozdziewski 2008), makes it 
possible to derive the analytic secular theory up to the order of 10 
in the semi-major axes ratio a. The accuracy of this model may be 
compared with the results of the numerical averaging of the pertur- 




Figure 17. Evolution of the mean orbits in the following system: mo = 
1 m , m\ = 100 mj, mi = 10 mj, a\ = 0.2 au, ai = 5 au, e\ = 0.45, 
0.001. coi = rc/2, co 2 = rc/2, iW = 67°.46, A = 0.215. Panels from the 
top to the bottom illustrate the evolution of e\ , <?2, *mut, Wi , CO2, respectively. 
The grey curves are for the classic secular model, the black curves are for 
the model including the relativistic corrections to the potential of the inner 
binary. 

bation ( Michtchenko & Malhotra 2004 ). Moreover, for a class of hi- 
erarchical systems considered in our paper, already the third-order 
model is precise enough to find and investigate the qualitative fea- 
tures of the system. In this work, we focus on three-body configura- 
tions, e.g., a star and two massive planets or a binary stars and one 
planet. "Fortunately", the second-order model is integrable, hence 
the octupole-level approximation might be considered as the first 
order perturbation to this analytically soluble case. This makes it 
possible to understand the sources of instabilities appearing in the 
full (non- averaged) model. 

The averaging over mean anomalies reduces the dynamics to 
a system having two DOF, which then may be investigated with 
the help of rich geometrical tools, like the Poincare cross section 
and the representative planes of initial conditions introduced in 
Michtchenko et al. (2006). Unfortunately, all additional corrections 
that increase the DOF number must be neglected here, as for in- 
stance the rotational quadrupole moment of the star and/or of the 
planets. This is the price that must be paid for the possibly global 
model of the dynamics. In the assumed range of a\ and masses, 
the relativistic "corrections" in fact compare with the Newtonian 
point-mass mutual interactions, and are much larger than other per- 
turbations, like the tidal and rotational distortions of the bodies. 
The analysis of the secular frequencies introduced by various cor- 
rections justify that the model only takes into account the PN per- 
turbation to the potential of star and the inner secondary. 

Having the two DOF model, we investigate the simplest class 
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Figure 13. Energy levels presented in the fP^-plane, and calculated for the following parameters a = 0.04,^ = 10, A = 0.215, io = 84°, a\ = 0.2 au,«2 = 
5.0 au,rao = 1 m©. A sequence of panels demonstrates the view of the phase space when masses of the secondary bodies decrease (their ratio is constant). 
From panel b) to panel d) (rai,ra2)[mj] are (100, 10), (18.5, 1.85), (10, 1). Shaded areas mark ranges of i mut larger than 60°, 75°, (the light- and dark- grey, 
respectively). 
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Figure 14. Dynamical maps for the classic, Newtonian model in the IPs-plane. The model parameters are a = 0.04, ju = 10,^t = 0.215, io = 84°, a\ 
0.2 au,#2 = 5.0 au,rao = 1 m©, m\ = 100 mj, = 10 mj. See the text and caption to Fig.|4]for the details. 
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Figure 15. The dynamical maps for the octupole model with relativistic corrections. Masses of the secondary bodies are mi = lOOmj, mi = 10mj. See the 
text and caption to Fig.|4]for more details. 
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Figure 16. The dynamical maps for the octupole model with relativistic corrections. Masses of the secondary bodies are m\ = 18.5mj, ra2 = 1.85mj. See the 
text and caption to Fig.|4]for more details. 



of solutions that are the equilibria. We focus on the Lidov-Kozai 
resonance (LKR) in systems characterized by large range of the 
mass ratio /u. This part extends the results derived for small jjl in 
( Migaszewski & Gozdziewski 2009a). We found that even much 
smaller outer body (planetary mass with respect to sub-stellar 
values of m\) moving in a wide, highly inclined (/ mut >^ 60°) or- 
bit may significantly perturb the inner orbit. In turn, the restricted 
model of the circumbinary planet, when we assume that the planet 



does not influence the binary, is not generally valid, even if the in- 
ner mass m\ is 100 times larger than the outer body m^. 

We also studied the parametric structure of families of partic- 
ular equilibria classified as IVa, Ilia, Illb, IVb+ in (Migaszewski 
|& Go zdziewski 2009a) for small a. Thanks to this assumption, the 
analytic model makes it possible to investigate the transition be- 
tween the planetary regime (small }S) and the circumbinary regime 
(ju ~ 50, 100). This study shows that solutions in the planetary prob- 
lem may disappear for large mass ratio. 
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A particularly interesting feature of the ocupole model is the 
appearance of the secular chaos. We found that if z mut exceeds the 
critical inclination / crit (^/,a), the long term evolution of the system 
may be strongly chaotic, leading to large amplification of the eccen- 
tricities. In the regular regions of the phase space, the mean angle 
coi librates around ±7t/2. In some parts of these regions also the 
second secular angle 002 librates around ±7t/2. The initial condi- 
tions satisfying / mut > / crit lead to strong amplification of the inner 
eccentricity e\. Simultaneously, for the same values of the angu- 
lar momentum of the system, we may observe strong amplifica- 
tion of £2 m some region of the phase space, with almost constant 
relative inclination of the orbits. This behaviour may be attributed 
to unstable equilibrium Illb emerging in the secular system. The 
amplification of e\ happens not only for librations of COi around 
±7t/2 (in the LK regime), but also and particularly when this an- 
gle behaves chaotically, varying in the whole range of [0,2tc]. The 
dynamical maps reveal that the primarily source of the chaotic mo- 
tions are the unstable equilibria and unstable periodic orbits in the 
full system, following the appearance of separatrices in the inte- 
grate, quadrupole-level model. 

Thanks to the simple analytic model, the influence of relativis- 
tic correction on the global secular dynamics of the problem 
might be clearly demonstrated. A simple proof of this influence is 
provided by the analysis of the equilibria in the perturbed model. 
The differences between the Newtonian and relativistic models are 
larger when the mutual interactions between the secondaries are 
weaker, e.g., when companion masses are smaller. Yet the dynam- 
ics are basically very simple up to the limit of the critical inclina- 
tion, when the first bifurcation of the origin {e\ = £2 = 0) occurs, 
and this feature of the dynamics is preserved in both models. 

We stress that although the analysis is done for specific, dis- 
crete mass ratios, the results are valid as far the assumptions of 
the averaging theorem are fulfilled and the corrections besides gen- 
eral relativity are negligible. We demonstrated that similarly to the 
co-planar problem, the global 3-D dynamics of the classic, Newto- 
nian model essentially depend only on the ratios of semi-major axes 
and masses of the secondaries. Hence, although we consider mostly 
the circumbinary configurations, the results may be are easily ex- 
trapolated to the "typical" planetary regime, investigated already 
( Michtchenko et al. 2006 ). Moreover, when the PN corrections are 
added to the model, the dynamics are much more complex. Still 
the global picture of the phase space is determined by the ratio of 
these corrections and the Newtonian mutual interactions. Hence, if 
the system scales down, while this ratio is roughly preserved, the 
structure of the phase space, determined by stationary solutions of 
the secular system should not essentially change. 

Our approach may be generalized for other perturbations, like 
the rotational and conservative tidal distortion of the bodies in the 
system. Unfortunately, in the most general case, the dimension of 
the hierarchical system cannot be generally reduced to two DOR 
Moreover, these perturbations lead to even more interesting and 
intriguing dynamics, which we investigated in the co-planar and 
spatial case (e.g., |Fabrycky & Tre maine 2007, Mardling 2007, 
|Ragozzine & Worf|2009||Migaszewski & Gozdziewski|2009bi rWe 
work on a global approach, suitable for the 3D systems, aiming to 
publish these results in future papers. 



ACKNOWLEDGEMENTS 

We thank the anonymous referee for comments that improved the 
manuscript. This work is supported by the Polish Ministry of Sci- 



ence and Higher Education, through grants N/N203/402739 and 
92/N-ASTROSIM/2008/0. CM is a recipient of the stipend of the 
Foundation for Polish Science (programme START, edition 2010). 



REFERENCES 

Adams F. C, Laughlin G., 2006, International Journal of Modern 
Physics D, 15,2133 

Arnold V. I., Kozlov V. V., Neishtadt A. I., 1993, Dynamical sys- 
tems III. Mathematical aspects of classical and celestial mechan- 
ics. Encyclopaedia of mathematical sciences, Springer Verlag 

Batygin K., Bodenheimer P., Laughlin G., 2009, ApJL, 704, L49 

Brouwer D., Clemence G. M., 1961, Methods of celestial mechan- 
ics. New York: Academic Press, 1961 

Brumberg V., 2007, Celestial Mechanics and Dynamical Astron- 
omy, 99, 245 

Eggenberger A., 2010, in K. Gozdziewski, A. Niedzielski, & 
J. Schneider ed., EAS Publications Series Vol. 42. pp 19-37 

Fabrycky D., Tremaine S., 2007, ApJ, 669, 1298 

Farago F, Laskar J., 2010, MNRAS, 401, 1189 

Ferraz-Mello S., ed. 2007, Canonical Perturbation Theories - De- 
generate Systems and Resonance Vol. 345 of Astrophysics and 
Space Science Library 

Ferrer S., Osacar C, 1994, Celestial Mechanics and Dynamical 
Astronomy, 60, 187 

Ford E. B., Kozinsky B., Rasio F. A., 2000, ApJ, 535, 385 

Gronchi G. F, Milani A., 1998, Celestial Mechanics and Dynam- 
ical Astronomy, 71, 109 

Harrington R. S., 1968, AJ, 73, 190 

KozaiY., 1962, AJ, 67, 579 

Krasinskii G. A., 1972, Celestial Mechanics, 6, 60 

Krasinskii G. A., 1974 Vol. 62 of IAU Symposium, pp 95-116 

Laskar J., 1990, Icarus, 88, 266 

Laskar J., 2000, Physical Review Letters, 84, 3240 

Lee M. H., Peale S. J., 2003, ApJ, 592, 1201 

Libert A.-S., Henrard J., 2007, Icarus, 191, 469 

Lidov M. L., 1962, Planetary and Space Science, 9, 719 

Lidov M. L., Ziglin S. L., 1976, Celestial Mechanics, 13, 471 

Mardling R. A., 2007, MNRAS, 382, 1768 

Mardling R. A., 2010, MNRAS, 407, 1048 

Michtchenko T. A., Ferraz-Mello S., Beauge C, 2006, Icarus, 
181,555 

Michtchenko T. A., Malhotra R., 2004, Icarus, 168, 237 
Migaszewski C, Gozdziewski K., 2008, MNRAS, 388, 789 
Migaszewski C, Gozdziewski K., 2009a, MNRAS, 395, 1777 
Migaszewski C, Gozdziewski K., 2009b, MNRAS, 392, 2 
Migaszewski C, Gozdziewski K., 2010, in K. Gozdziewski, 
A. Niedzielski, & J. Schneider ed., EAS Publications Series 
Vol. 42, .pp 385-391 
Murray C. D., Dermott S. F, 2000, Solar System Dynamics. Cam- 
bridge Univ. Press 
Rabl G., Dvorak R., 1988, A&A, 191, 385 
Ragozzine D., Wolf A. S., 2009, ApJ, 698, 1778 
Richardson D. L., Kelly T. J., 1988, Celestial Mechanics, 43, 193 
Takeda G., Kita R., Rasio F. A., 2009, in IAU Symposium Vol. 253 

of IAU Symposium, . pp 181-187 
Tamuz, O., et al., 2008, A&A, 480, L33 

Sidlichovsky M., Nesvorny D., 1996, Celestial Mechanics and 
Dynamical Astronomy, 65, 137 



© 2008 RAS, MNRAS 000,[l][T8] 



